subroutine AxisyPossion2DT31
    use Mesh2DT3_m
    use FES2DT31_m
    use SMatrix_m
    use IterSSolver_m
    implicit none
    type(Mesh2DT3) :: mesh2D
    type(FES2DT31) :: fs
    type(SMatrix) :: HeatStiff
    type(IterSSolver) :: solver
    real*8, allocatable :: x(:), RHS(:)
    real*8, allocatable :: xExact(:), error0(:)
    real*8 :: L2Norm, L2RError
    integer :: BDLabel(2) = (/2, 4/)

    call mesh2D.square(100, 10, 2)
    call mesh2D.scale(1.0D0, 0.1D0)
    call mesh2D.translate(1.0D0, 0.0D0)

    call fs.initFES(mesh2D)
    call HeatStiff.initSMatrix(fs%nDof)
    allocate(x(fs%nDof), RHS(fs%nDof), xExact(fs%nDof), error0(fs%nDof))
    x = 0.0D0
    RHS = 0.0D0
    
    call fs.assembleHeatStiff(1.0D0, AXISYMMETRIC, HeatStiff)
    call fs.assembleHeatSource(10.0D0, AXISYMMETRIC, RHS)
    
    call fs.applyBCZero(HeatStiff, RHS, BDLabel)
    call solver.initIterSSolver()
    call solver.PCGSSOR(HeatStiff, RHS, x)

    call fs.valueFromfunc(uExact, xExact)

    error0 = x - xExact
    L2Norm = fs.computeL2Norm(xExact)
    L2RError = fs.computeL2Norm(error0) / L2Norm
    write(*, *) "L2 RError", L2RError

contains

real*8 function uExact(xy)
    implicit none
    real*8, intent(in) :: xy(2)
    real*8 :: x, y
    x = xy(1); y = xy(2)
    uExact = -2.5D0 * x**2 + 7.5D0 * log(real(x)) / log(2.0D0) + 2.5D0
end function



end subroutine
